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Confined  Two-Phase  Incompressible  Flows 

T.  S.  Mautner 
NCCOSC  RDTE  Division* 

San  Diego,  CA  92152-6171 

Abstract  translating  container.  The  MAC  method^'^  has  been 

extended  to  allow  for  variable  fluid  properties  and 
The  variable  density,  incompressible  Navier-Stokes  surface  tension.  Results  of  the  flow  field  and  the 

equations  are  solved  using  a  technique  based  upon  the  interface  motion  will  be  presented  along  with  estimates 

Marker-and-Cell  (MAC)  method.  The  flow  equations  of  the  conservation  of  mass  and  stability  requirements, 

have  been  supplemented  with  the  level  set  approach  for 

interface  tracking,  and  the  required  modifications  to  the  Nomenclature 

MAC  method  are  presented.  This  formulation  has  been 

used  to  solve  several  example  problems,  for  closed  D  divergence 

containers,  involving  single  and  two-phase  flows,  Fr  Froude  number  =  U/Lg 

various  boundary  conditions  and  several  types  of  wall  F^  force  due  to  surface  tension 

motion.  g  acceleration  due  to  gravity 

h  grid  cell  dimension  -  (Ax,Ay) 

Introduction  L  characteristic  length 

n  unit  normal  vector 

A  survey  of  recent  literature  indicates  a  strong  p  pressure 

interest  in  the  computation  of  incompressible  flows  Re  Reynolds  number  =  UL/v 

with  density  variations,  free  surfaces  and  bubbles  rising  t  time 

and  interacting  with  one  another  and  with  walls’’*"'®'’^'’^.  At  time  step 

An  integral  part  of  these  problems  is  the  computation  of  At^  time  step  due  to  convection 

the  interface  between  the  various  fluid  regions.  In  Af,  time  step  due  to  viscous  terms 

contrast  to  Lagrangian  methods  and  the  use  of  mass-less  u  fluid  velocity  (u,v) 

particles^  various  interface  tracking  and  volume-of-  u„  magnitude  for  translating  box 

fluid  methods  have  been  developed  and  provide  an  u,^  wall  velocity 

effective  way  to  compute  multi-phase  flows  (for  U  characteristic  velocity 

examples,  see  Refs.  4-7,  9-11,13-16).  We  Weber  number  =  pU2L/o 

x,y  Cartesian  coordinates 

Recently  Osher  and  Sethian"  introduced  a  6  delta  function 

"Hamilton-Jacobi"  level  set  formulation  of  the  p  fluid  density 

equations  for  propagating  interfaces.  The  level  set  p,  density  ratio  =  pj/p, 

approach  implicitly  captures  the  interface  location,  as  v  kinematic  viscosity 

the  zero  level  of  a  smooth  function,  rather  than  p  dynamic  viscosity 

explicitly  tracking  the  front.  With  this  method  steep  p,  viscosity  ratio  =  p^/p, 

fluid  gradients  can  be  handled  as  well  as  the  merging  x  viscous  stress  tensor 

and  breaking  of  fluid  regions.  The  level  set  approach  0^  rotation  angle 

has  been  used  to  study  fluid  instabilities  in  compressible  o  surface  tension  coefficient 

flows  and  incompressible  flows  with  flame  burning,  k  curvature 

shear  layers,  driven  cavities  and  bubbles.  w  frequency 

cOj,  over-relaxation  parameter 

This  paper  applies  the  level  set  approach"-'^  to  a  interface  thickness 

problems  of  fluids  in  closed  containers  with  various  (()  level  set  function 

forms  of  wall  motion  such  as  the  driven  cavity  and  level  set  function  at  previous  time 

*This  paper  is  declared  work  of  the  U.S.  Government  and 
is  not  subject  to  copyright  protection  in  the  United  States. 
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Governing  Equations 

The  equations  of  motion  are  the  variable 
density,  incompressible  Navier-Stokes  equations' 


V-u=0  (1) 


(pu),  +  V-(puu)  = 
-Vp  +  Vr  +  pg  +  Fj 


V-T  =  pV^u  +  Vp'^  (6) 


In  general  the  velocities  computed  from  Eqn 
(2)  will  not  satisfy  the  continuity  equation  Eqn.  (1).  The 
SOLA  method  does  not  solve  an  elliptic  equation  for 
the  pressure;  however,  the  divergence  free  condition  is 
obtained  by  adjusting  the  cell  pressures.  The  pressures 
are  determined  iteratively,  and  the  divergence  D  in  each 
cell  is  computed  using  the  most  current  values  of  the 
velocity  components.  The  pressure  change,  Ap, 
required  to  make  the  divergence  equal  zero  is 


For  immiscible  liquids,  the  density  and  viscosity  are 
constant  along  particle  paths,  thus  they  can  be 
computed  using 


Ap 


-co„  D 

_ 0 _ 

2  At  {(Ax)-'^  +  (Ay)-2} 


(7) 


pj  +  (u-^p  =  0  (3) 


p,  +  (ti-^p  =  0  (4) 


Using  characteristic  length,  velocity,  density  and 
viscosity  parameters,  the  non-dimensional  form  of  the 
momentum  equation  is 


(pit),  =  -V(puu)  - 

Fr 

+  (-Vp  +  -iv-T  +  Fj  ) 


(5) 


The  Navier-Stokes  equations  are  solved  using 
a  modified  form  of  the  SOLA^  formulation  of  the  MAC 
method  of  Harlow  and  Welch’.  The  method  was 
modified  to  account  for  variable  fluid  properties,  and 
the  finite-difference  expressions  for  the  equations 
utilized  the  staggered  MAC  grid  presented  in  Fig.  1. 
Even  though  this  formulation  is  for  constant  density 
fluids,  the  pressure  and  density  are  considered  separate 
terms  instead  of  a  single  term,  i.e.  the  kinematic 
pressure  (p/p),  as  in  the  majority  of  incompressible 
models.  Additionally,  the  viscosity  is  considered  a 
variable  in  the  viscous  term  formulation  given  by 


and  the  new  pressure,  at  the  k(th)  iteration,  is 


k.l  k  * 

Pio  =  Pij  ^  ^P 

(8) 

The  velocity  components,  at  the  cell  faces,  are 
adjusted  to  reflect  the  change  in  pressure  using  the 
following  relationships 

k*l  k  At  A 

Uii  =  Uii  +  — Ap 

•j  .d 

(9) 

k*l  k  At  A 

^i-ij  =  “i-ij  - 

(10) 

k*l  k  At  A 

%  -  -'ij  * 

(11) 

k*l  k  At  A 

''iA-.  ■ 

(12) 

where  the  density  is  evaluated  at  the  cell  faces  using,  for 
example. 

,  _  (Pi.i  ^  Pi*i,i) 

n.i/2j  2 

(13) 
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Convergence  of  the  pressure  solution  can  be  accelerated 
by  multiplying  Ap  by  an  over-relaxation  parameter,  a)„, 
where 

Since  the  density  and  viscosity  change  sharply 
across  an  interface,  conventional  finite-difference 
techniques  may  not  be  able  to  handle  the  steep 
derivatives  that  would  occur  at  these  points.  Thus,  to 
capture  the  fluid  interface,  the  level  set  technique' 
will  be  used.  The  level  set  function  (({))  will  be  zero 
along  the  interface,  positive  in  one  fluid  and  negative  in 
the  other.  The  level  set  function  is  initialized  as  the 
signed  normal  distance  to  the  interface  and  has  the  form 


(j),  +  (u-tO(j)  =0  (14) 


This  equation  will  transport  the  zero  level  set  of  ({) 
exactly  as  the  fluid  interface  moves.  Also,  since  (J)  is  a 
smooth  function,  it  is  more  easily  solved  than  Eqns.  (3) 


and  (4)  for  the  density  and  viscosity.  Now,  in  this 

formulation,  the  density  and  viscosity  are  defined  by 

[l 

if 

4)>o 

P  =  ' 

P2/P1 

if 

4x0  (15) 

[(Pi+PjVP: 

if 

4)=o 

’  1 

if 

4)>o 

P  =  j 

P2/P1 

if 

4><0  (16) 

(P,+P2)/P2 

if 

4)=o 

The  level  set  technique  may  develop  jumps 
when  interfaces  interact,  thus  the  level  set  function  will 
be  reinitialized  to  keep  the  distance  function  a  distance 
function.  First,  to  prevent  instabilities  at  the  interface, 
especially  for  large  density/viscosity  ratios,  the  density 
and  viscosity  can  be  smoothed  according  to 

p  =  (pj+p2)/2p,  p  =  (p,+P2)/2p; 

Ap  =  .(p,-p2)/2p;  Ap  =  (p,-p2)/2p, 


f  1  if  (t»0 

P(4))=  Pz/Pi  (17) 

[  p  +  Ap s i n (7i(t)/2a)  otherwise 


and 


f  1  if  (j)>0 

p(4>)=]P2(Pi  'f  4><0  (18) 

[p+Ap  sin(::4>/2a)  otherwise 


It  should  be  noted  that  both  the  density  and  viscosity 
are  smoothed  across  the  interface  region  in  contrast 
to  Sussman  et.  al'^  who  did  not  smooth  the  viscosity 
in  this  manner. 

Reinitialization  of  the  distance  function  is 
performed  using  the  equation 

4).  =  S((t)„)(l-/^^}  (19) 


where 


(20) 


Applying  the  technique  given  by  Sussman  et. 
al.'^,  the  following  expressions  are  used  to  solve  Eqn. 
(19) 


a=(4)ij-<|)i.,j)/Ax  c=((t)i.-(j)i..,)/Ay 


if  (j)=0 


G((j))  = 


1 -(  max(b,0)^+min(a,0)^  if  <t)<0 
+max(d,0)^+min(c,0)^  )'^ 
l-(  max(a,0)^+min(b,0)^  if  (j)>0 
+max(c,0)^+min(d,0)^  )'^ 


(22) 
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Using  these  expressions,  Eqn  (19)  now  has  the  form 
(1)>=+At  S((l)„)  G((j))  (23) 

and  an  iterative  technique  is  used  to  drive  the  solution 
to  steady  state. 

The  surface  tension  term  in  Eqn  (2)  is 


In  the  translating  container  case,  an  additional  term  was 
added  to  the  x  momentum  equation  to  account  for  the 
coordinate  translation.  The  term  is 

-pu^cosin(cL)t)  (29) 


and  for  container  rotation,  the  gravity  term  becomes 


Fj  =  0  K  6(d)  h  (24) 


per„,sin(a)t)-l- 

F 


(30) 


and  can  be  expressed  in  terms  of  the  distance  function 
as 


1 

We 


K((J))  6(4))  V4) 


(25) 


where  the  curvature  is 


k(4))=V-|-^  . 
[lV4)l 


(26) 


and  the  delta  function  6(4))  can  be  approximated  by 


Finally,  at  the  solid  walls,  the  level  set  function 
maintained  a  zero  normal  gradient.  Flowever,  other 
boundary  conditions  could  have  been  applied. 

It  is  well  known  that  the  explicit  formulations 
have  time  step  limitations  based  upon  the  CFL 
condition  for  the  convective  terms,  the  viscous  terms 
and  the  stiff  source  terms  due  to  gravity  and  surface 
tension.  The  criteria  used  herein  are  based  on  the 
recommendations  in  Refs.  3,  6,  and  13.  For  example, 
one  can  use 


At 


min /Ax  Ay| 
Ci  .  lul  ’  Ivl  j 


(31) 


and 


6(4))  = 


2 

— {l  +cos(7t4)/a)}/a 

I4)l<a 

2 

0 

u  (27) 

other¬ 

,,  min  1  .  pReh-^ 

At,  s  ——(const  -t- - ) 

Q  1  V  J 

wise 

^  A 

(32) 


In  smoothing  the  fluid  properties  and  in  the  surface 
tension  term,  a  reasonable  value  for  the  interface 
thickness  is  a52Ax.  The  examples  in  this  paper  were 
obtained  using  a=0  and  no  surface  tension  effects  were 
included. 

The  boundary  conditions  applied  to  the 
examples  in  this  paper  are  both  the  no-slip  and  slip 
conditions  at  the  solid  walls.  For  the  driven  cavity 
problem,  the  top  wall  has  a  prescribed  velocity  of 


then 

At  i  min  (At„  ,  AtJ  (33) 

where  the  constant  was  suggested  to  be  (3/14).'^ 
However,  a  value  of  approximately  0.043-  0.13  was 
used  to  obtain  the  time  steps  for  the  example  cases.  The 
time  step  data,  obatained  using  Eqn.  (32)  exclusively, 
are  shown  in  Table  1 . 


Uw3„sin(cot)  (28) 
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Finally,  the  solution  algorithm  can  be 
summarized  as  follows; 

Step  1.  Initialize  (j)  to  be  the  signed  normal 
distance  to  the  interface. 

Step  2.  Initialize  the  fluid  regions  with  the 
appropriate  density  and  viscosity  values. 

Step  3.  Compute  the  (u,v)  velocity  components 
using  Eqns.  (5)-(6). 

Step  4.  Compute  the  change  in  pressure  to 
obtain  a  divergence  free  flow  and  update  the 
(u,v)  velocity  components,  Eqns.  (7)-(12). 

Step  5.  Construct  a  new  distance  function  by 
solving  Eqn.  (14)  to  steady  state. 

Step  6.  Adjust  the  density  and  viscosity  in  the 
fluid  regions,  Eqns.  (17)-(18). 

Repeat  steps  3-6  for  each  At. 

Example  Calculations 

The  computational  method  described  above 
was  applied  to  several  problems  involving  confined 
flows.  The  results  for  a  driven  cavity,  translating 
container,  rotating  container  and  a  density  driven  flow 
will  be  given  below.  The  input  parameters  are 
summarized  in  Table  1,  and  the  arrangement  of 
variables  in  a  typical  computational  cell  is  shown  in 
Figure  1.  It  should  be  noted  that,  in  all  computations, 
surface  tension  was  not  imposed,  g/(FR)^=l,  and  the 
interface  between  the  two  fluid  densities  is  located  at 
the  middle  of  the  container  at  t=0. 

1.  Driven  Cavity 

The  constant  fluid  property  formulation  of 
the  method  was  tested  using  a  driven  cavity  having 
from  5x5  to  200x200  rectangular  grids.  The  results 
compared  very  favorably  with,  for  example,  SOLA’ 
results.  Figure  2  presents  the  velocity  and  pressure 
fields  for  a  cavity  with  Re=1000  with  a  constant  lid 
velocity,  and  the  results  agree  well  with  computations 
made  using  a  vorticity-stream  function  method. 

The  next  case  considered  was  the  driven 
cavity  having  a  sinusoidal  lid  velocity  with  frequency 
co=0.63  and  Re=100.  The  frequency  of  oscillation 
produces  two  complete  cycles  in  approximately  t=20 
non-dimensional  time  units.  The  results  for  the  constant 
fluid  property  case  are  a  given  in  Figs.  3  and  4  where 
the  no-slip  and  slip  boundary  conditions  were  used. 


Both  the  velocity  field  and  pressure  contour  data  show 
the  expected  shift  of  the  fluid's  center  of  rotation  due  to 
the  time-dependent  lid  velocity.  At  the  points  of 
maximum  lid  velocity,  the  results  compare  well  with  the 
uniform  lid  velocity  case.  Comparison  of  the  data  in 
Figs  3  and  4  shows  the  shift  of  the  velocity’s  center  of 
rotation  with  the  change  in  boundary  conditions. 

Next,  using  a  small  density  ratio,  p,=0.9, 
constant  viscosity  and  co=0.315,  the  effect  of  grid 
density  was  examined.  The  results  for  the  50x50  grid 
are  presented  in  Fig.  5  and  the  100x100  grid  results  are 
given  in  Fig.  6.  It  is  clear  from  the  time  histories  of  the 
interface  location  and  the  velocity  field  that  the  solution 
has  not  converged.  Higher  grid  densities  should  be  used 
to  determine  solution  convergence.  Another 
comparison,  using  the  100x100  grid,  was  made  between 
the  results  obtained  using  the  no-slip  and  slip  boundary 
conditions.  The  results  found  in  Figs.  6  and  7  once 
again  show  the  differences  in  the  computed  results  for 
the  no-slip  and  slip  boundary  conditions.  These 
differences  are  not  unexpected  and  should  remind  one 
of  the  importance  of  using  appropriate  boundary 
conditions.  The  results  also  show  that  the  level  set 
formulation  has  no  problem  following  the  interface  as 
it  moves  along  the  cavity  wall  and  when  it  wraps  up  due 
to  the  rotating  velocity  field. 

The  final  driven  cavity  example  uses  density 
and  viscosity  ratios  of  0.5  and  0.2  respectively  and  a 
frequency  a)=0.63.  The  computed  results  are  given  in 
Fig.  8  and  show  that,  with  this  higher  density  ratio,  the 
velocity  field  has  a  large  rotation  confined  to  the  lower 
density  region  (top)  while  a  smaller  velocity  rotation 
occurs  in  the  heavy  fluid  region  (bottom).  In  this  case 
the  interface  had  only  a  small  oscillatory  motion.  When 
one  increases  the  density  ratio  to  10,  the  fluid  motion  is 
almost  entirely  confined  to  the  lighter  (top)  fluid  region 
with  little  or  no  interface  motion.  The  p,=0.1  results  are 
not  shown  here. 

2.  Translating  Container 

The  next  case  considered  was  that  of  a 
container,  having  a  2:1  aspect  ratio,  translating, 
sinusoidally  in  the  horizontal  direction,  with  frequency 
(0=0.63  and  magnitudes  Uo=1.0  and  0.25,  see  Eqn.  (29). 
Both  cases  used  the  same  density  (pj=0.5)  and  viscosity 
(p,=0.2)  ratios,  and  the  time  histories  of  the  interface 
motion  and  velocity  fields  are  shown  in  Figs.  9  and  10. 
From  the  results  of  the  higher  magnitude  of  oscillation 


5 


case  shown  in  Fig.  9  (u„=1.0),  the  ability  of  the  level 
set  method  to  follow  highly  convoluted  interface  motion 
is  clearly  demonstrated.  The  interface  contours  look 
somewhat  like  that  of  a  breaking  wave.  In  contrast,  the 
Uo=0.25  case  shows  a  gentle  sloshing  of  the  fluid. 

3.  Rotating  Container 

Rotation  of  a  rectangular  container  about  its 
mid-point,  aspect  ratio  of  4:1,  was  considered  using  the 
same  density  ratios  as  above  but  now  considering  a 
maximum  rotation  angle  of  30  degrees  with  frequencies 
of  rotation  (0=0.315  and  (o=0.628.  The  results  of  the 
interface  motion  and  velocity  field,  over  time,  are 
plotted  in  Figs.  11  and  12.  As  before,  the  level  set 
method  follows  the  interface  at  all  times  and  indicates 
a  wave  breaking  like  shape  for  the  co=0.315  case.  The 
data  in  Fig.  12,  co=0.628,  indicates  that  the  flow  has 
assumed  a  gentle  sloshing  pattern. 

4.  Density  Driven  Flow 

The  final  example  considered  was  a 
Rayleigh-Taylor  instability  problem.  A  density  ratio  of 
p,=i).5  and  a  viscosity  ratio  of  iir=0.2  were  used  where 
the  heavy  fluid  was  located  in  the  upper  half  of  the 
domain.  The  problem  was  started  with  the  application 
of  a  small  negative  velocity  (v=-0.25)  along  the  middle 
1/3  of  the  interface.  The  time  evolution  of  the  interface 
is  given  in  Fig.  13  and  shows  the  typical  behavior  found 
by  many  researchers.  The  slight  asymmetric  result  is 
due  to  the  uneven  number  of  grid  cells  on  each  side  of 
the  mid-point  along  the  interface.  As  in  the  many  recent 
bubble  computations,  the  interface  moves  smoothly 
along  the  wall  at  large  times. 

5.  Interface  Resolution 

To  illustrate  the  transition  of  fluid 
properties  across  the  interface,  selected  portions  of  the 
interface  from  the  rotating  container  at  t=17  (Fig.  1 1) 
and  the  Rayleigh-Taylor  problem  at  t=10  (Fig.  13)  are 
shown  in  Fig.  14.  It  can  be  seen  that  for  the  “zero” 
interface  thickness  used  in  these  examples,  the  density 
transitions  over  one  grid  cell  with  the  zero  of  the  level 
set  function  falling  within.  The  interface  locations  and 
these  density  contours  were  obtained  using  the  data 
computed  at  the  cell  centers  shown  in  Fig.  1.  If  one 
plots  interface  and  density  contours  using  node  point 
averaged  data,  the  interface  will  be  smeared  over 
approximately  3  grid  cells. 


Stability  and  Conservation 

Use  of  the  stability  criteria  given  in  Eqns. 
(31)-  (33)  provided  reasonable  values  for  the  time  step 
in  this  explicit,  finite-difference  computation  provided 
that  the  constant  value  in  Eqn.  (32)  was  reduced  from 
3/14‘^  to  approximately  0.043-0.13.  The  time  step 
values  given  in  Table  1  are  conservative;  however,  the 
computation  remained  stable  over  long  computation 
times.  The  computations  also  showed  that  the  density 
and  viscosity  ratios  play  an  important  role  in 
determining  the  appropriate  time  step.  Also,  an  estimate 
of  mass  conservation  (area  in  two  dimensions)  was 
made  by  simply  counting  the  number  of  cell  center 
points  in  the  positive  and  negative  regions  of  the  level 
set  function  (]).  The  maximum  error  results  are  given  in 
Table  1  and  were  determined  from  the  maximum  area 
ratios  greater  than  and  less  than  the  nominal  area  ratio 
at  t=0.  It  is  encouraging  that,  even  for  the  “zero” 
thickness  interface,  the  maximum  errors  are  less  than 
10%  with  the  largest  errors  occurring  during  large 
distortions  of  the  interface. 

Conclusion 

The  computed  results  presented  above 
indicate  successful  modification  of  the  incompressible 
MAC/SOLA  method  to  include  variable  density  flows 
augmented  with  the  level  set  approach  for  interface 
tracking.  The  finite-difference  expressions  followed  the 
work  of  Daly*  and  Hirt  et.  al?  for  the  Navier-Stokes 
equations  and  a  second  order  formulation  for  the  level 
set  method.  The  example  problems  demonstrated  the 
ability  of  the  method  to  track  interface  motion  using  a 
“zero”  thickness  interface  (1  grid  cell  width)  for  density 
ratios  up  to  10.  Higher  density  ratios  could  also  be 
handled  using  some  smoothing  of  the  density  and 
viscosity  across  the  interface,  for  example  Eqns.  (17)- 
(18).  The  effects  of  grid  density,  wall  motion  and 
boundary  conditions  on  the  solution  were  also 
presented.  Potential  uses  of  this  method  could  be  the 
study  of  breaking  waves  and  the  transport  of  immiscible 
material. 
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Table  1 .  Input  parameters  and  error  estimates  for  example  problems. 


Example 

Fig. 

■  No. 

Grid  Size 
h 

■ 

F, 

B.C.* 

At 

^wall 

%  Max 
Error’ 

Driven 

Cavity 

2 

200x200 

0.005 

1.0 

1.0 

No-Slip 

0.001 

G)=0 

1 .0 

Driven 

Cavity 

3 

50x50 

0.02 

1.0 

1.0 

No-Slip 

0.005 

w=0.630 

Uw.n=l-0 

Driven 

Cavity 

■ 

50x50 

0.02 

1.0 

1.0 

Slip 

0.005 

05=0.630 

Uw»,=1.0 

Driven 

Cavity 

5 

50x50 

0.02 

0.9 

1.0 

No-Slip 

0.005 

(0=0.315 

Uv«ii=l-0 

6.62/2.30 

Driven 

Cavity 

6 

100x100 

0.01 

0.9 

1.0 

No-Slip 

0.001 

(0=0.315 
11,^=  1.0 

7.96/0.19 

Driven 

Cavity 

■ 

100x100 

0.01 

0.9 

1.0 

Slip 

0.001 

(0=0.315 

u,^=1.0 

2.55/2.71 

Driven 

Cavity 

8 

50x50 

0.02 

0.5 

0.2 

No-Slip 

0.005 

(0=0.630 

u,^=1.0 

1.54/0.86 

Translating 

Container 

9 

100x50 

0.02 

0.5 

0.2 

No-Slip 

0.005 

(0=0.63 

u„=1.0 

3.45/1.63 

Translating 

Container 

10 

100x50 

0.02 

0.5 

0.2 

No-Slip 

0.005 

(0=0.63 

u„=0.25 

0.96/0.77 

Rotating 

Container 

11 

200x50 

0.02 

0.5 

0.2 

No-Slip 

0.005 

(0=0.315 

0^=30” 

7.58/2.50 

Rotating 

Container 

12 

200x50 

0.02 

0.5 

0.2 

No-Slip 

0.005 

(i)=0.628 

0ro,=3O‘’ 

0.77/0.38 

Rayleigh- 

Taylor 

13 

100x200 

0.01 

0.5 

0.2 

Slip 

0.001 

2.84/5.16 

Notes:  1.  B.C.  -  Boundary  Condition 

2.  Frequency  used  in:  a)  Eqn(28)  for  Figs.  2-8,  b)  Eqn(29)  for  Figs.  9-10; 
c)  Eqn.(30)  for  Figs.  1 1-12. 

3.  Max  error  refers  to  %  greater/less  than  nominal  area  ratio  for  interface. 
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Figure  2.  Velocity  and  pressure  fields  at  t=4.9  for  a  driven  cavity,  Re=1000,  no¬ 
slip  boundary  conditions  and  constant  fluid  properties. 
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L  0.040 
K  0.0M 
J  0.032 
1  0.021 
H  0.024 
Q  0.020 
f  0.016 
E  0.012 
0  0.001 
C  0.004 
8  0.000 
A  -0.004 
9  -0.006 

1  -0.012 

7  -0.016 

6  -0.020 
S  -0.024 
4  -0.028 

3  -0.032 

2  -0.036 

1  -0.040 


L  0.040 
K  0.036 
J  0.032 

1  0.028 

H  0.024 
Q  0.020 
F  0.016 
E  0.012 
D  0.006 
C  0.004 
B  0.000 
A  -0.004 
9  -0.006 

8  -0.012 
7  -0.016 

6  -0.020 
S  -0.024 
4  -0.026 

3  -0.032 

2  -0.036 

1  -0.040 


Figure  3.  Velocity  and  pressure  fields  at  t=5  and  t=10  for  a  driven  cavity,  Re=100, 
no-slip  boundary  conditions  and  constant  fluid  properties. 


L  0.100 
K  0.090 
J  0.060 
I  0.070 
H  0.060 

a  0.050 

F  0.040 
E  0.030 
0  0.020 
C  0.010 
B  0.000 
A  -0.010 
9  -0.020 

•  -0.030 

7  -0.040 

6  -0.060 
9  -0.000 

4  -0.070 

3  *0.060 

2  -0.090 

1  -0.100 


Figure  4.  Velocity  and  pressure  fields  at  t=5  for  a  driven  cavity,  Re=100,  slip 
boundary  conditions  and  constant  fluid  properties. 
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Figure  5.  Interface  location  and  velocity  field  for  a  driven  cavity,  Re=100,  no-slip 
boundary  conditions,  p  =  0.9,  Pj=l,0  and  grid  size  of  50x50. 


Figure  6.  Interface  location  and  velocity  field  for  the  driven  cavity,  Re=100, 
no-slip  boundary  conditions,  Pr=0.9,  Mr=1.0  and  grid  size  of  100x100. 
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Figure  8.  Interface  location  and  velocity  field  for  the  driven  cavity,  Re=100, 
no-slip  boundary  conditions,  Pr=0.5,  Pr=0.2  and  grid  size  of  50x50. 


Figure  10.  Interface  location  and 
Uo=0.25,  pr=0.5  and  \i=0.2. 


field  for  a  translating  box,  Re=100 


Figure  11.  Interface  location  and  velocity  filed  for  a  rotating  box,  Re=100, 
0„,=3O*sin(O.315-t),  p  =0.5  and  p,=0.2. 


Figure  12.  Interface  location  and  velocity  filed  for  a  rotating  box,  Re=100, 
0^=3O*sin(O.628-t),  p  =0.5  and  ^,=0-2. 


fluid  Rayleigh-Taylor  instability, 
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Figure  14.  Examples  of  density  interface  contours  for  the  rotating  container  and 
the  Rayieigh-Taylor  examples. 
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